#load libraries 
library(pacman)
p_load("tidyverse","ggpubr", "ggspatial", "cowplot", "here","gridExtra","rgeos","rgdal","rmapshaper", "MASS", "nlme", "mgcv","data.table", "kableExtra","gtsummary","Hmisc","arsenal","imputeTS","zoo") 
source(here::here("COVID","print_table_flextable.R"))

Introduction

This file includes all R code for the analysis carried out for the paper “Impact of COVID-19 pandemic on mental health medical service utilization in Australian youth: an interrupted time series analysis”.

Item table

Items<- read_csv(here::here("Data", "Item code.csv")) %>% 
    mutate(`Group Description`=str_to_sentence(`Group Description`)) %>% 
    dplyr::select(`MBS Group`,`Group Description`,`MBS Requested Item`) %>% 
    group_by(`MBS Group`,`Group Description`) %>% 
    mutate(Items = paste0(`MBS Requested Item`, collapse = ", ")) %>% 
    dplyr::select( -`MBS Requested Item`) %>% 
    unique() %>% 
    mutate(`Group Description`=str_replace(`Group Description`,"Gp", "GP")) %>% 
    mutate(`Group Description`=str_replace(`Group Description`,"gp", "GP")) %>% 
    mutate(`Group Description`=str_replace(`Group Description`,"Covid", "COVID"))


print_table(Items, align = c("l", "l", "l"), cap = "MBS items included", capname = "Items", colwidth = c(1,5,10))
(#tab:Items)MBS items included

Import data

The data include service uptake count data by SA3, year and population subgroups. Other avaliable inforamtion include SA3 region type (regional or metropolitan), state (Victoria vs. other states), estimated residential population (ERP) and the 2016 Socio-economic Indexes for Areas (SEIFA) scores for each SA3 (estimated based on population weighted average of SA2 scores) were used in this analysis. The measure of SEIFA used was the Index of Relative Socio-economic Advantage and Disadvantage (IRSAD).

dta <- readRDS( here("COVID","COVID_cleaned.rds"))
map <- readRDS(here("Data", "Shapefile_SA3.rds"))

MBS service rates by SA3

This heat map displays the annual number of services accessed per 1000 people for the combined youth population in 200 for each SA3 in Australia. Separate figures show the heat maps specifically for Greater Brisbane, Greater Sydney and Greater Melbourne.

library(raster)
# state boundary
State <- raster::aggregate(map, "STE_NAME16")

# calculate annual service access rates
rates <- dta %>%
  filter( Population == "All young people") %>%
  group_by(SA3_CODE16, Year) %>%
  summarise(Services = sum(Services), ERP = mean(ERP)) %>%
  mutate(Rates = Services * 1000 / ERP) %>%
  mutate(SA3_CODE16 = as.character(SA3_CODE16))
map@data <- map@data %>%
  left_join(rates[rates$Year == 2020, ], by = "SA3_CODE16")

Aus <- ggplot() +
  # plot pologon
  ggspatial::annotation_spatial(map) +
  layer_spatial(map, aes(fill = Rates), colour = "gray68") +
  layer_spatial(State, colour = "gray20", fill = NA) +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(
    location = "bl", which_north = "true",
    pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = c(0.15, 0.85),
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.border = element_blank(),
  ) +
  labs(fill = "Number of services \n per 1000 population") +
  coord_sf(xlim = c(109, 165), ylim = c(-9, -45))

Mel <- ggplot() +
  # plot pologon
  ggspatial::annotation_spatial(map) +
  layer_spatial(map, aes(fill = Rates), colour = "gray68") +
  scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  coord_sf(xlim = c(144.45, 145.55), ylim = c(-37.5, -38.3))

Syd <- ggplot() +
  # plot pologon
  ggspatial::annotation_spatial(map) +
  layer_spatial(map, aes(fill = Rates), colour = "gray68") +
  scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  coord_sf(xlim = c(150.67, 151.32), ylim = c(-33.65, -34.15))


Bri <- ggplot() +
  # plot pologon
  ggspatial::annotation_spatial(map) +
  layer_spatial(map, aes(fill = Rates), colour = "gray68") +
  scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  coord_sf(xlim = c(152.65, 153.5), ylim = c(-27.05, -27.75))

arrows <- data.frame(xs = c(0.62, 0.71, 0.74), xe = c(0.79, 0.79, 0.79), ys = c(0.23, 0.33, 0.5), ye = c(0.18, 0.43, 0.67))

ggdraw(Aus) +
  # add AUD map
  draw_plot(Bri, x = 0.75, y = 0.55, width = 0.25, height = 0.25) +
  draw_plot(Syd, x = 0.75, y = 0.3, width = 0.25, height = 0.25) +
  draw_plot(Mel, x = 0.75, y = 0.05, width = 0.25, height = 0.25) +
  geom_segment(aes(x = xs, y = ys, xend = xe, yend = ye),
    data = arrows,
    arrow = arrow(), lineend = "round"
  )

Imputation

Series were removed from the dataset if they had fewer than 15 non-missing data points. The remaining missing values were imputed using linear interpolation, generating integers between 0 and 10.

# extract missing
dta <- dta %>%
  group_by(SA3_CODE16, Population) %>%
  mutate(missing = sum(is.na(Services))) %>% 
  ungroup()

# check missing
# table(dta$missing)

set.seed(12345)
dta_imputed<-dta %>%
  #keep if the series had less than 15 missing data points
  filter( missing <= 15) %>% 
  # Remove few data points for newly established areas where ERP was too small 
  filter(ERP>5)%>% 
  # fill missing with random number between 0 and 9
  mutate(Services=na_random(Services,lower_bound= -0.5, upper_bound=9.5)) %>% 
  mutate(Services=round(Services)) 
  
#prepare data for GLM 
dta_imputed <- dta_imputed %>% 
  mutate(YearQuarter=as.numeric(YearQuarter),
         COVID = as.numeric(Year==2020),
         ID = paste0(SA3_CODE16,"_",Population)) %>% 
   arrange(ID,YearQuarter)

# checking
# hist(dta$Services[dta$Services<50])
# hist(dta_imputed$Services[dta_imputed$Services<50])
# length(unique(dta_imputed$SA3_CODE16))

Line plot

The line plot below presents the quarterly number of services delivered per 1000 people for each SA3 by population subgroup, along with the mean rate at each quarter for that subpopulation.

dta %>%
  mutate(YearQuarter = as.Date(YearQuarter)) %>%
  group_by(YearQuarter,Population) %>% 
  mutate(Mean = mean(Rates,na.rm = TRUE) ) %>% 
  ungroup() %>% 
  ggplot(aes(x = YearQuarter, y = Rates, group = SA3_CODE16)) +
  geom_line(alpha = 0.5, col = "grey60", size=0.2) +
  geom_line( aes(x=YearQuarter, y=Mean),  col = "red") +
  facet_wrap(vars(Population), ncol = 1) +
  scale_x_date(date_labels = "%Y", date_breaks = "1 year") +
  labs(x = "Year", y = "Number of services per 1000 population") +
  theme_bw() +
  theme(
    strip.background=element_rect(fill="white")
  )

IRSAD distribution

IRSAD<-dta_imputed %>% 
  dplyr::select(SA3_NAME,IRSAD) %>% 
  distinct()
summary(IRSAD)
##            SA3_NAME       IRSAD       
##  Adelaide City :  1   Min.   : 748.6  
##  Adelaide Hills:  1   1st Qu.: 943.3  
##  Albany        :  1   Median : 983.1  
##  Albury        :  1   Mean   : 995.5  
##  Alice Springs :  1   3rd Qu.:1044.8  
##  Armadale      :  1   Max.   :1166.6  
##  (Other)       :325
a<-IRSAD %>% 
  ggplot(aes(IRSAD)) +
  stat_ecdf(geom = "step") + 
  theme_bw()  + 
  scale_y_continuous(labels=scales::percent)+
  scale_x_continuous(limits = c(740,1170)) +
  labs(y="Percentage")

b<-IRSAD %>% 
  mutate(IRSAD=plyr::round_any(IRSAD,30)) %>% 
  ggplot(aes(IRSAD)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) + 
  theme_bw()  + 
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(limits = c(740,1170)) +
  labs(y="Percentage")

ggpubr::ggarrange(a, b,
          labels = c("A", "B"),nrow = 2)

# Descriptive table

This table describes the median (IQR) number of services delivered per 1000 people at an SA3 level for 2019 and 2020 for each subgroup, as well as the percentage change in these values.

tbl_dta<-dta %>% 
  filter(Year>=2019) %>% 
  group_by(SA3_CODE16,Population,Year) %>% 
  summarise(Services=sum(Services), ERP=mean(ERP)) %>% 
  mutate(Rates=Services*1000/ERP) %>%
  ungroup() %>% 
  dplyr::select(Population,Year,Rates,SA3_CODE16) %>% 
  pivot_wider(names_from=Year,values_from=Rates) %>% 
  mutate(`Percentage increase`= ( `2020`- `2019`)*100/ `2019`)


tbl_dta %>% 
  dplyr::select(-SA3_CODE16) %>% 
  tbl_summary(
             by = Population,
             missing = "no",
             digits = `Percentage increase` ~ 1,
          ) 
Characteristic Age 12-17, N = 3321 Age 18-25, N = 3321 Females, N = 3321 Males, N = 3321 All young people, N = 3321
2019 572 (428, 688) 696 (552, 850) 837 (661, 1,008) 438 (346, 543) 639 (496, 760)
2020 622 (470, 787) 793 (618, 997) 984 (779, 1,231) 443 (360, 569) 717 (551, 889)
Percentage increase 12.6 (4.0, 21.0) 15.7 (9.9, 21.8) 20.4 (14.0, 28.0) 2.8 (-2.9, 8.5) 14.0 (8.4, 20.3)

1 Median (IQR)

Boxplot

These boxplots show the distribution of the number of services accessed per 1000 people for each quarter between 2019 and 2020 by population subgroup and year.

dta %>% 
  filter(Year>=2019) %>% 
  mutate(Quarter=as.factor(paste0("Q", Quarter)),
         Year=as.factor(Year)) %>% 
  mutate(Rates=Services*1000/ERP) %>%
  ggplot(aes(x=Quarter, y=Rates, col=Year)) +
  geom_boxplot() +theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "bottom",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#0072B2","#E69F00")) +
  labs(y="Rates per 1000 population")

dta %>% 
  filter(Year>=2019) %>% 
  mutate(Quarter=as.factor(paste0("Q", Quarter)),
         Year=as.factor(Year)) %>% 
  mutate(Rates=Services*1000/ERP) %>%
  dplyr::select(Quarter,Rates,Year,Population,SA3_NAME) %>% 
  pivot_wider(names_from="Population",values_from = "Rates") %>% 
  dplyr::select(-SA3_NAME) %>% 
  tbl_strata(
    strata = Quarter,
    .tbl_fun =
      ~ .x %>%
        tbl_summary(by = Year, missing = "no")
  )
Characteristic Q1 Q2 Q3 Q4
2019, N = 3321 2020, N = 3321 2019, N = 3321 2020, N = 3321 2019, N = 3321 2020, N = 3321 2019, N = 3321 2020, N = 3321
Age 12-17 132 (98, 161) 134 (100, 168) 148 (111, 181) 148 (111, 192) 152 (114, 185) 173 (132, 226) 134 (99, 165) 163 (119, 211)
Age 18-25 167 (133, 206) 176 (137, 219) 173 (136, 215) 197 (152, 245) 183 (142, 223) 219 (169, 272) 169 (131, 204) 202 (158, 250)
Females 199 (155, 244) 209 (163, 260) 214 (168, 260) 243 (187, 306) 223 (173, 264) 281 (216, 343) 198 (155, 239) 260 (203, 321)
Males 108 (81, 132) 106 (82, 131) 111 (88, 138) 110 (90, 140) 118 (91, 145) 122 (97, 156) 109 (85, 131) 108 (85, 140)
All young people 151 (117, 184) 155 (122, 195) 161 (126, 198) 174 (135, 219) 168 (131, 203) 200 (155, 249) 153 (119, 184) 180 (142, 230)

1 Median (IQR)

State comparision

Conversely, these boxplots depict the distribution of the yearly number of services delivered per 1000 people for Victoria and the rest of Australia by subpopulation and year. The yearly rates are obtained similarly to the formula for the heat maps above.

dta %>% 
  filter(Year>=2019 ) %>% 
  mutate(Quarter=as.factor(paste0("Q", Quarter)),
         Year=as.factor(Year)) %>% 
  group_by(SA3_CODE16,Population,Year,State) %>% 
  summarise(Services=sum(Services), ERP=mean(ERP)) %>% 
  mutate(Rates=Services*1000/ERP) %>%
  ungroup() %>% 
  ggplot(aes(x=State, y=Rates, col=Year)) +
  geom_boxplot() +theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "bottom",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#0072B2","#E69F00")) +
  labs(y="Rates per 1000 population")

dta %>% 
  filter(Year>=2019 ) %>% 
  mutate(Quarter=as.factor(paste0("Q", Quarter)),
         Year=as.factor(Year)) %>% 
  group_by(SA3_CODE16,Population,Year,State) %>% 
  summarise(Services=sum(Services), ERP=mean(ERP)) %>% 
  mutate(Rates=Services*1000/ERP) %>%
  ungroup() %>% 
  dplyr::select(State,Rates,Year,Population,SA3_CODE16) %>% 
  pivot_wider(names_from="Year",values_from = "Rates") %>%
  mutate(`Percentage increase`= (`2020`-`2019`)*100/`2019`) %>% 
  dplyr::select(-SA3_CODE16) %>% 
  tbl_strata(
    strata = Population,
    .tbl_fun =
      ~ .x %>%
        tbl_summary(by = State, missing = "no",
                    digits = `Percentage increase` ~ 1)
  )
Characteristic Age 12-17 Age 18-25 Females Males All young people
Other, N = 2661 Victoria, N = 661 Other, N = 2661 Victoria, N = 661 Other, N = 2661 Victoria, N = 661 Other, N = 2661 Victoria, N = 661 Other, N = 2661 Victoria, N = 661
2019 553 (413, 680) 636 (508, 733) 682 (525, 839) 765 (610, 870) 813 (639, 999) 921 (751, 1,089) 435 (335, 536) 474 (405, 550) 612 (482, 757) 695 (560, 817)
2020 619 (457, 783) 655 (511, 855) 761 (595, 972) 886 (739, 1,071) 963 (744, 1,205) 1,155 (933, 1,405) 439 (344, 556) 465 (408, 596) 685 (540, 874) 796 (659, 984)
Percentage increase 13.3 (5.3, 21.3) 8.9 (2.1, 17.4) 14.3 (8.8, 20.2) 20.9 (15.3, 26.2) 19.4 (13.4, 27.1) 25.6 (16.3, 31.1) 2.9 (-3.2, 8.6) 2.5 (-2.7, 7.7) 13.1 (8.2, 19.8) 16.9 (10.1, 22.3)

1 Median (IQR)

Using data from 2015-2020

Indiviudal model

The interrupted time series models for the mean number of services at time \(t\), \(\mu_t\), took the form:

\[\log(\mu_t) = \beta_0 + \beta_1t + \beta_2COVID_t + \beta_{3}Quarter_{t} + \log(ERP_t)\]

,where \(COVID_t\) is the interruption jump coefficient, \(Quarter_{t}\) is the quarterly seasonal effect term and \(ERP_t\) is the estimated resident population at time \(t\).

#initialize results
results<-list()
count=1

# loop indiviudal time sereis
for(id in unique(sort(dta_imputed$ID))) {
    
    temp <-dta_imputed %>% 
      filter(ID==id & Year >= 2015)
    #check for evidence of residual autocorrelation 
    model <- glm(Services ~ YearQuarter  + COVID + Quarter +  offset(log(ERP)) ,
                 data = temp, family = "quasipoisson")
    res<-resid(model, type='pearson')
    test<-Box.test(res, lag=10, fitdf=0)
    
    
    # when p is large will not control for autocorrelation 
    if(test$p.value>0.05){
      #convergence issues with the default optimiser in a few series
      if(tryCatch(model<-MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),  
                             random = ~ 1|SA3_NAME, 
                             data =temp, family = "quasipoisson"), 
              error = function(c) "error")=="error")
        {
        model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
              random = ~ 1|SA3_NAME, temp, family = "quasipoisson",
              control = lmeControl(opt = "optim"))
      }
    } else {
      count=count+1
      #convergence issues with the default optimiser in a few series
      if(tryCatch(model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
                             random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
                             correlation = corAR1()), error = function(c) "error")=="error") {
            model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
              random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
              control = lmeControl(opt = "optim"),
              correlation = corAR1())
            }
    }
  
  #save results 
  result <- coef(summary(model)) 
  result <-cbind(rownames(result),as_tibble(result))
  result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
  result$Population <-unique(temp[temp$ID==id,]$Population)
  results[[id]]<-result 
}


#extract SA3 information 
sa3 <- dta %>% 
    dplyr::select(SA3_CODE16,SA3_NAME,IRSAD,Region, State)  %>% 
    unique()


results<-as_tibble(rbindlist(results,idcol="ID")) %>% 
    mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>% 
    left_join(sa3,by="SA3_CODE16")

names(results)[2]<-"Terms"

#fraction of models with autocorrelation issues 
count/length(unique(sort(dta_imputed$ID)))
## [1] 0.1212675
write.csv(results %>% filter(Terms=="COVID"),"results.csv")

Plots

This scatter plot shows the estimated relative risk against the IRSAD score for each subgroup. Superimposed is locally weighted fitted smoothed line and the 95% CI for the prediction at that IRSAD score.

COVID<-results %>% 
  filter(Terms=="COVID") %>% 
  mutate(RR=exp(Value)) 

COVID %>% 
  #3 sa3 excluded at the tail
  dplyr::filter( IRSAD > 810) %>% 
  ggplot(aes(x=IRSAD, y=RR)) +
  geom_point( alpha=0.6, col="#0072B2", size=1)+
  geom_smooth(col="red",span=1) +
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "bottom",
        strip.background=element_rect(fill="white"))

These boxplots compare the distribution of the estimated relative risks by state (Victoria vs. rest of Australia) and region (metropolitan vs. reginal) for each subpopulation.

state<-COVID %>% 
  ggplot(aes(x=State, y=RR, colour=State)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#56B4E9","#0072B2"))

region<- COVID %>% 
  ggplot(aes(x=Region, y=RR, colour=Region)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#E69F00", "#D55E00"))

ggpubr::ggarrange(state,region, 
          labels = c("A", "B"), nrow = 2)

These two scatter plots compare the estimated relative risks within each population variable (age group and sex) and their relationships with the colouring variable, IRSAD score. The dashed line represents equality between the estimates for each of the sexes/age-groups.

library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Males" | Population=="Females") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=Females>Males)  %>%
  ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for males", y="RR for females") +
  theme(legend.position="none") 

b<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Age 12-17" | Population=="Age 18-25") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=`Age 18-25` >`Age 12-17`)  %>%
  ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for age 12-17", y="RR for age 18-25")+
  theme(legend.position="none") 

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

legend<-a +
      theme(legend.position="bottom") +
      theme(legend.key.height= unit(0.5, 'cm'),
        legend.key.width= unit(1.2, 'cm'))
      
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
             legend, nrow=2,heights=c(10, 1))

Pooled effect

This forest plot portrays the average estimated relative risk (95% CI) for each subgroup as estimated by a random effects meta-analysis model.

library(metafor)
extract_meta<-function(data){
  a<-rma(Value, Std.Error, data=data,method = "REML", test = "knha")
  ret<-exp(c(a$b,a$ci.lb, a$ci.ub ))
  names(ret)<-c( "RR", "CIL","CIH")
  ret
}

res <- results %>% 
   filter(Terms=="COVID") %>% 
   group_by(Population) %>%
   group_modify(~ broom::tidy(extract_meta(data = .x))) %>% 
   pivot_wider(names_from=names,values_from=x) %>% 
   mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
                                    format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3),")")) 
res$y_loc=5:1



#main plot in the middle 
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
        geom_point() +
        geom_errorbarh(aes(xmin = CIL, xmax = CIH),  height = 0.3) +
        ylab(NULL) +
        scale_y_continuous(limits = c(0.8,5.5)) +
        ggtitle(" ") +
        xlab("Estimated RR (95% CI)") +
        scale_colour_manual(values = colours ) +
        theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
              axis.text.y = element_blank(),
              axis.ticks.y = element_blank(),
              panel.border = element_blank(),
              axis.line.x = element_line(colour = "black"),
              panel.background = element_blank(),
              panel.grid = element_blank(),
              plot.margin = margin(0, 0, 0, 0, "cm"),
              legend.position = "none")  + 
        geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50") 

#base plot for adding text 
tab_base <- ggplot(res, aes(y = y_loc)) +
        ylab(NULL) + xlab("") + 
        scale_y_continuous(limits = c(0.8,5.5)) +
        theme(plot.title = element_text( size = 12), ## centering title on text
              axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
              axis.line = element_blank(),
              axis.text.y = element_blank(),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              plot.margin = margin(0, 0, 0.3, 0, "cm"),
              legend.position = "none",
              panel.background = element_blank(),
              panel.grid = element_blank()) +
         scale_colour_manual(values = colours )  +
        labs(x="  ")

# Predictors
tab1 <- tab_base + 
        geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
        scale_x_continuous(limits = c(0,0.4)) +
        ggtitle("Outcomes") +
        theme(plot.title = element_text(hjust = 0.3)) 

# OR
OR1 <- tab_base +
        geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) + 
        ggtitle(" OR (95% CI)")+
        theme(plot.title = element_text(hjust = 0.5))


gridExtra::grid.arrange(tab1, plot1, OR1,
             layout_matrix = matrix(c(rep(1,2),
                                      rep(2,6),
                                      rep(3,3)), nrow = 1))

Meta regression

Random effects meta-regression was used to compile the interruption effect estimates and identify their relationships with the predictors IRSAD score, state (rest of Australia vs. Victoria) and area type (metropolitan vs. regional). The meta-regression model took the form:

\[\hat{\beta}^{i}_2 = \alpha_0 + \alpha_1Regional_i+\alpha_2IRSAD_i+ \alpha_3Victoria_i + u_{i} + \epsilon_{i}\]

, where \(\hat{\beta}^{ij}_2\) is the estimated coefficient for the COVID-19 impact for SA3 \(i\) in a specific population. \(Regional_i\) is the binary area type variable (baseline metropolitan), \(IRSAD_i\) is the IRSAD score for SA3 \(i\) centred on the mean IRSAD score for Australia and expressed in 100 units, \(Victoria_i\) is the binary state variable (baseline rest of Australia), \(\epsilon_{i}\) is the sampling error and \(u_{i}\) is the error due to heterogeneity between SA3s in true interruption effect size.

library(metafor)
reg_meta<-function(data){
  a<-rma(yi=Value, sei=Std.Error, data=data,method = "REML", 
         mods = ~  IRSAD + State + Region, test = "knha") 
  browser()
  print( unique(data$population))
  print(summary(a))
  ret<-data.frame(cbind(exp(cbind(a$b, a$ci.lb,a$ci.ub)), a$pval))
  names(ret)<-c( "RR", "CIL","CIH", "p")
  ret$Variable<-c("Intercept","IRSAD","Victoria", "Regional" )
  return(ret)
}

res <- results %>% 
   filter(Terms=="COVID") %>% 
   mutate(IRSAD= (IRSAD-995.5)/100) %>% 
   mutate(population=Population)%>%
   group_by(Population) %>%
   group_modify(~reg_meta(data = .x)) %>% 
   mutate('RR ' = format(round(RR, digits=3), nsmall=3),
          "95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3)),
          'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
   dplyr::select(-RR,-CIL,-CIH,-p) %>% 
   pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>% 
   pivot_wider(names_from = Population,values_from=Value)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Age 12-17
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  250.9328  -501.8655  -491.8655  -472.9773  -491.6762   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0071 (SE = 0.0009)
## tau (square root of estimated tau^2 value):             0.0840
## I^2 (residual heterogeneity / unaccounted variability): 67.04%
## H^2 (unaccounted variability / sampling variability):   3.03
## R^2 (amount of heterogeneity accounted for):            20.88%
## 
## Test for Residual Heterogeneity:
## QE(df = 323) = 971.5387, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 19.5891, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt           0.0337  0.0090   3.7440  0.0002   0.0160  0.0515  *** 
## IRSAD             0.0746  0.0105   7.1266  <.0001   0.0540  0.0952  *** 
## StateVictoria    -0.0240  0.0145  -1.6615  0.0976  -0.0525  0.0044    . 
## RegionRegional    0.0223  0.0143   1.5534  0.1213  -0.0059  0.0505      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Age 18-25
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  354.9734  -709.9469  -699.9469  -681.0432  -699.7582   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0030 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0550
## I^2 (residual heterogeneity / unaccounted variability): 59.08%
## H^2 (unaccounted variability / sampling variability):   2.44
## R^2 (amount of heterogeneity accounted for):            44.50%
## 
## Test for Residual Heterogeneity:
## QE(df = 324) = 800.0453, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 41.8378, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt           0.0374  0.0064  5.8261  <.0001  0.0248  0.0500  *** 
## IRSAD             0.0495  0.0074  6.6672  <.0001  0.0349  0.0641  *** 
## StateVictoria     0.0962  0.0109  8.8435  <.0001  0.0748  0.1176  *** 
## RegionRegional    0.0304  0.0103  2.9648  0.0033  0.0102  0.0506   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Females
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  346.9755  -693.9511  -683.9511  -665.0628  -683.7618   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0029 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0542
## I^2 (residual heterogeneity / unaccounted variability): 54.42%
## H^2 (unaccounted variability / sampling variability):   2.19
## R^2 (amount of heterogeneity accounted for):            42.61%
## 
## Test for Residual Heterogeneity:
## QE(df = 323) = 716.2943, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 37.7070, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval   ci.lb   ci.ub 
## intrcpt           0.0833  0.0065  12.7335  <.0001  0.0705  0.0962  *** 
## IRSAD             0.0664  0.0077   8.6598  <.0001  0.0513  0.0815  *** 
## StateVictoria     0.0636  0.0111   5.7289  <.0001  0.0418  0.0854  *** 
## RegionRegional    0.0339  0.0105   3.2279  0.0014  0.0132  0.0546   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Males
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  293.9581  -587.9161  -577.9161  -559.0124  -577.7274   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0044 (SE = 0.0006)
## tau (square root of estimated tau^2 value):             0.0667
## I^2 (residual heterogeneity / unaccounted variability): 63.89%
## H^2 (unaccounted variability / sampling variability):   2.77
## R^2 (amount of heterogeneity accounted for):            11.88%
## 
## Test for Residual Heterogeneity:
## QE(df = 324) = 892.8698, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 9.7570, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval    ci.lb    ci.ub 
## intrcpt          -0.0494  0.0078  -6.3100  <.0001  -0.0648  -0.0340  *** 
## IRSAD             0.0449  0.0089   5.0314  <.0001   0.0273   0.0624  *** 
## StateVictoria     0.0141  0.0125   1.1269  0.2606  -0.0105   0.0386      
## RegionRegional    0.0165  0.0124   1.3292  0.1847  -0.0079   0.0410      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] All young people
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 331; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  382.8617  -765.7234  -755.7234  -736.7736  -755.5364   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0026 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0515
## I^2 (residual heterogeneity / unaccounted variability): 60.75%
## H^2 (unaccounted variability / sampling variability):   2.55
## R^2 (amount of heterogeneity accounted for):            39.45%
## 
## Test for Residual Heterogeneity:
## QE(df = 327) = 835.6417, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 327) = 35.2839, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt           0.0373  0.0059  6.2989  <.0001  0.0256  0.0489  *** 
## IRSAD             0.0591  0.0069  8.6269  <.0001  0.0457  0.0726  *** 
## StateVictoria     0.0498  0.0100  4.9931  <.0001  0.0302  0.0694  *** 
## RegionRegional    0.0266  0.0095  2.8057  0.0053  0.0079  0.0452   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
res$Variable[c(2:3,5:6,8:9,11:12)]<-""


print_table(res, align = c("l", "l",rep("c",5)),  colwidth = c(2,3,rep(4.4,5)))
## <caption> (\#tab:)</caption>

Sensitivity analysis using data from 2011

All results below are generated similarly as before for longer time series from 2011.

Individual model

#initialize results
results<-list()
count=1

# loop indiviudal time sereis
for(id in unique(sort(dta_imputed$ID))) {
    
    temp <-dta_imputed %>% 
      filter(ID==id )
    #check for evidence of residual autocorrelation 
    model <- glm(Services ~ YearQuarter  + COVID + Quarter +  offset(log(ERP)) ,
                 data = temp, family = "quasipoisson")
    res<-resid(model, type='pearson')
    test<-Box.test(res, lag=10, fitdf=0)
    
    
    # when p is large will not control for autocorrelation 
    if(test$p.value>0.05){
      #convergence issues with the default optimiser in a few series
      if(tryCatch(model<-MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),  
                             random = ~ 1|SA3_NAME, 
                             data =temp, family = "quasipoisson"), 
              error = function(c) "error")=="error")
        {
        model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
              random = ~ 1|SA3_NAME, temp, family = "quasipoisson",
              control = lmeControl(opt = "optim"))
      }
    } else {
      count=count+1
      #convergence issues with the default optimiser in a few series
      if(tryCatch(model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
                             random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
                             correlation = corAR1()), error = function(c) "error")=="error") {
            model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)), 
              random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
              control = lmeControl(opt = "optim"),
              correlation = corAR1())
            }
    }
  
  #save results 
  result <- coef(summary(model)) 
  result <-cbind(rownames(result),as_tibble(result))
  result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
  result$Population <-unique(temp[temp$ID==id,]$Population)
  results[[id]]<-result 
}


#extract SA3 information 
sa3 <- dta %>% 
    dplyr::select(SA3_CODE16,SA3_NAME,IRSAD,Region, State)  %>% 
    unique()


results<-as_tibble(rbindlist(results,idcol="ID")) %>% 
    mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>% 
    left_join(sa3,by="SA3_CODE16")

names(results)[2]<-"Terms"

#fraction of models with autocorrelation issues 
count/length(unique(sort(dta_imputed$ID)))
## [1] 0.4582572

Plots

COVID<-results %>% 
  filter(Terms=="COVID") %>% 
  mutate(RR=exp(Value)) 

COVID %>% 
  #3 sa3 excluded at the tail
  dplyr::filter( IRSAD > 810  ) %>% 
  ggplot(aes(x=IRSAD, y=RR)) +
  geom_point( alpha=0.6, col="#0072B2", size=1)+
  geom_smooth(col="red",span=1) +
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "bottom",
        strip.background=element_rect(fill="white"))

state<-COVID %>% 
  ggplot(aes(x=State, y=RR, colour=State)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#56B4E9","#0072B2"))

region<- COVID %>% 
  ggplot(aes(x=Region, y=RR, colour=Region)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#E69F00", "#D55E00"))

ggpubr::ggarrange(state,region, 
          labels = c("A", "B"), nrow = 2)

library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Males" | Population=="Females") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=Females>Males)  %>%
  ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for males", y="RR for females") +
  theme(legend.position="none") 

b<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Age 12-17" | Population=="Age 18-25") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=`Age 18-25` >`Age 12-17`)  %>%
  ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for age 12-17", y="RR for age 18-25")+
  theme(legend.position="none") 

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

legend<-a +
      theme(legend.position="bottom") +
      theme(legend.key.height= unit(0.5, 'cm'),
        legend.key.width= unit(1.2, 'cm'))
      
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
             legend, nrow=2,heights=c(10, 1))

Pooled effect

res <- results %>% 
   filter(Terms=="COVID") %>% 
   group_by(Population) %>%
   group_modify(~ broom::tidy(extract_meta(data = .x))) %>% 
   pivot_wider(names_from=names,values_from=x) %>% 
   mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
                                    format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3),")")) 
res$y_loc=5:1

#main plot in the middle 
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
        geom_point() +
        geom_errorbarh(aes(xmin = CIL, xmax = CIH),  height = 0.3) +
        ylab(NULL) +
        scale_y_continuous(limits = c(0.8,5.5)) +
        ggtitle(" ") +
        xlab("Estimated RR (95% CI)") +
        scale_colour_manual(values = colours ) +
        theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
              axis.text.y = element_blank(),
              axis.ticks.y = element_blank(),
              panel.border = element_blank(),
              axis.line.x = element_line(colour = "black"),
              panel.background = element_blank(),
              panel.grid = element_blank(),
              plot.margin = margin(0, 0, 0, 0, "cm"),
              legend.position = "none")  + 
        geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50") 

#base plot for adding text 
tab_base <- ggplot(res, aes(y = y_loc)) +
        ylab(NULL) + xlab("") + 
        scale_y_continuous(limits = c(0.8,5.5)) +
        theme(plot.title = element_text( size = 12), ## centering title on text
              axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
              axis.line = element_blank(),
              axis.text.y = element_blank(),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              plot.margin = margin(0, 0, 0.3, 0, "cm"),
              legend.position = "none",
              panel.background = element_blank(),
              panel.grid = element_blank()) +
         scale_colour_manual(values = colours )  +
        labs(x="  ")

# Predictors
tab1 <- tab_base + 
        geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
        scale_x_continuous(limits = c(0,0.4)) +
        ggtitle("Outcomes") +
        theme(plot.title = element_text(hjust = 0.3)) 

# OR
OR1 <- tab_base +
        geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) + 
        ggtitle(" OR (95% CI)")+
        theme(plot.title = element_text(hjust = 0.5))


gridExtra::grid.arrange(tab1, plot1, OR1,
             layout_matrix = matrix(c(rep(1,2),
                                      rep(2,6),
                                      rep(3,3)), nrow = 1))

Meta regression


res <- results %>% 
   filter(Terms=="COVID") %>% 
   mutate(IRSAD= (IRSAD-995.5)/100) %>% 
   group_by(Population) %>%
   group_modify(~reg_meta(data = .x)) %>% 
   mutate('RR ' = format(round(RR, digits=3), nsmall=3),
          "95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3)),
          'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
   dplyr::select(-RR,-CIL,-CIH,-p) %>% 
   pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>% 
   pivot_wider(names_from = Population,values_from=Value)

res$Variable[c(2:3,5:6,8:9,11:12)]<-""


print_table(res, align = c("l", "l",rep("c",5)),  colwidth = c(2,3,rep(4.4,5)))

Sensitivity analysis exclude autocorrelation

All results below are generated similarly as before for interrupted time series model without AR(1) structure estimated using the glm function.

Individual model

results<-list()

for(id in unique(sort(dta_imputed$ID))) {
  temp <-dta_imputed %>% 
      filter(ID==id & Year >= 2015)
  model <- glm(Services ~ Year  + COVID + Quarter +  offset(log(ERP)) , data = temp, family = "quasipoisson")
  result <- coef(summary(model)) 
  result <-cbind(rownames(result),as_tibble(result))
  result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
  result$Population <-unique(temp[temp$ID==id,]$Population)
  results[[id]]<-result 
}


#combine results
results<-as_tibble(rbindlist(results,idcol="ID")) %>% 
    mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>% 
    left_join(sa3,by="SA3_CODE16")

names(results)[2:4]<-c("Terms", "Value","Std.Error")

Plots

COVID<-results %>% 
  filter(Terms=="COVID") %>% 
  mutate(RR=exp(Value)) 

COVID %>% 
  #3 sa3 excluded at the tail
  dplyr::filter( IRSAD > 810  ) %>% 
  ggplot(aes(x=IRSAD, y=RR)) +
  geom_point( alpha=0.6, col="#0072B2", size=1)+
  geom_smooth(col="red",span=1) +
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "bottom",
        strip.background=element_rect(fill="white"))

state<-COVID %>% 
  ggplot(aes(x=State, y=RR, colour=State)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#56B4E9","#0072B2"))

region<- COVID %>% 
  ggplot(aes(x=Region, y=RR, colour=Region)) +
  geom_boxplot()+
  facet_wrap( ~ Population,ncol = 5)+ 
  theme_bw()  +
  facet_wrap( ~ Population,ncol = 5)+  
  theme(legend.position = "none",
        strip.background=element_rect(fill="white"))+
  scale_color_manual(values=c("#E69F00", "#D55E00"))

ggpubr::ggarrange(state,region, 
          labels = c("A", "B"), nrow = 2)

library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Males" | Population=="Females") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=Females>Males)  %>%
  ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for males", y="RR for females") +
  theme(legend.position="none") 

b<-COVID %>% 
  dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
  filter(Population=="Age 12-17" | Population=="Age 18-25") %>% 
  pivot_wider(names_from=Population,values_from=RR) %>% 
  mutate(compare=`Age 18-25` >`Age 12-17`)  %>%
  ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
  geom_jitter( alpha=0.8,  size=2)+
  theme_bw()  +
  scale_color_gradientn(colours = pal)  +
  scale_x_continuous(limits = c(0.6, 1.4)) +
  scale_y_continuous(limits = c(0.6, 1.4)) +
  geom_abline(slope = 1, col="gray50",linetype = 2)+ 
  labs(x="RR for age 12-17", y="RR for age 18-25")+
  theme(legend.position="none") 

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

legend<-a +
      theme(legend.position="bottom") +
      theme(legend.key.height= unit(0.5, 'cm'),
        legend.key.width= unit(1.2, 'cm'))
      
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
             legend, nrow=2,heights=c(10, 1))

Pooled effect

res <- results %>% 
   filter(Terms=="COVID") %>% 
   group_by(Population) %>%
   group_modify(~ broom::tidy(extract_meta(data = .x))) %>% 
   pivot_wider(names_from=names,values_from=x) %>% 
   mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
                                    format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3),")")) 
res$y_loc=5:1



#main plot in the middle 
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
        geom_point() +
        geom_errorbarh(aes(xmin = CIL, xmax = CIH),  height = 0.3) +
        ylab(NULL) +
        scale_y_continuous(limits = c(0.8,5.5)) +
        ggtitle(" ") +
        xlab("Estimated RR (95% CI)") +
        scale_colour_manual(values = colours ) +
        theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
              axis.text.y = element_blank(),
              axis.ticks.y = element_blank(),
              panel.border = element_blank(),
              axis.line.x = element_line(colour = "black"),
              panel.background = element_blank(),
              panel.grid = element_blank(),
              plot.margin = margin(0, 0, 0, 0, "cm"),
              legend.position = "none")  + 
        geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50") 

#base plot for adding text 
tab_base <- ggplot(res, aes(y = y_loc)) +
        ylab(NULL) + xlab("") + 
        scale_y_continuous(limits = c(0.8,5.5)) +
        theme(plot.title = element_text( size = 12), ## centering title on text
              axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
              axis.line = element_blank(),
              axis.text.y = element_blank(),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              plot.margin = margin(0, 0, 0.3, 0, "cm"),
              legend.position = "none",
              panel.background = element_blank(),
              panel.grid = element_blank()) +
         scale_colour_manual(values = colours )  +
        labs(x="  ")

# Predictors
tab1 <- tab_base + 
        geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
        scale_x_continuous(limits = c(0,0.4)) +
        ggtitle("Outcomes") +
        theme(plot.title = element_text(hjust = 0.3)) 

# OR
OR1 <- tab_base +
        geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) + 
        ggtitle(" OR (95% CI)")+
        theme(plot.title = element_text(hjust = 0.5))


gridExtra::grid.arrange(tab1, plot1, OR1,
             layout_matrix = matrix(c(rep(1,2),
                                      rep(2,6),
                                      rep(3,3)), nrow = 1))

Meta regression

res <- results %>% 
   filter(Terms=="COVID") %>% 
   mutate(IRSAD= (IRSAD-995.5)/100) %>% 
   group_by(Population) %>%
   group_modify(~reg_meta(data = .x)) %>% 
   mutate('RR ' = format(round(RR, digits=3), nsmall=3),
          "95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
                                    format(round(CIH, digits = 3),nsmall=3)),
          'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
   dplyr::select(-RR,-CIL,-CIH,-p) %>% 
   pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>% 
   pivot_wider(names_from = Population,values_from=Value)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  247.4643  -494.9285  -484.9285  -466.0403  -484.7393   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0074 (SE = 0.0009)
## tau (square root of estimated tau^2 value):             0.0861
## I^2 (residual heterogeneity / unaccounted variability): 68.96%
## H^2 (unaccounted variability / sampling variability):   3.22
## R^2 (amount of heterogeneity accounted for):            22.36%
## 
## Test for Residual Heterogeneity:
## QE(df = 323) = 1020.5763, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 21.3740, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt           0.0316  0.0092   3.4470  0.0006   0.0136  0.0497  *** 
## IRSAD             0.0786  0.0106   7.4014  <.0001   0.0577  0.0995  *** 
## StateVictoria    -0.0254  0.0147  -1.7243  0.0856  -0.0543  0.0036    . 
## RegionRegional    0.0220  0.0146   1.5132  0.1312  -0.0066  0.0507      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  348.7922  -697.5843  -687.5843  -668.6806  -687.3957   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0032 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0569
## I^2 (residual heterogeneity / unaccounted variability): 61.52%
## H^2 (unaccounted variability / sampling variability):   2.60
## R^2 (amount of heterogeneity accounted for):            43.50%
## 
## Test for Residual Heterogeneity:
## QE(df = 324) = 847.5023, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 40.4836, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt           0.0378  0.0066  5.7557  <.0001  0.0249  0.0508  *** 
## IRSAD             0.0503  0.0076  6.6052  <.0001  0.0353  0.0653  *** 
## StateVictoria     0.0966  0.0112  8.6598  <.0001  0.0747  0.1185  *** 
## RegionRegional    0.0315  0.0105  3.0002  0.0029  0.0109  0.0522   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  349.4359  -698.8718  -688.8718  -669.9836  -688.6826   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0030 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0545
## I^2 (residual heterogeneity / unaccounted variability): 55.25%
## H^2 (unaccounted variability / sampling variability):   2.23
## R^2 (amount of heterogeneity accounted for):            42.91%
## 
## Test for Residual Heterogeneity:
## QE(df = 323) = 732.5255, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 38.8790, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval   ci.lb   ci.ub 
## intrcpt           0.0838  0.0065  12.8567  <.0001  0.0710  0.0966  *** 
## IRSAD             0.0665  0.0076   8.6920  <.0001  0.0514  0.0815  *** 
## StateVictoria     0.0655  0.0111   5.8995  <.0001  0.0437  0.0874  *** 
## RegionRegional    0.0323  0.0105   3.0947  0.0021  0.0118  0.0529   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  294.2438  -588.4875  -578.4875  -559.5838  -578.2988   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0046 (SE = 0.0006)
## tau (square root of estimated tau^2 value):             0.0682
## I^2 (residual heterogeneity / unaccounted variability): 65.46%
## H^2 (unaccounted variability / sampling variability):   2.89
## R^2 (amount of heterogeneity accounted for):            12.17%
## 
## Test for Residual Heterogeneity:
## QE(df = 324) = 923.1324, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 10.1745, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     tval    pval    ci.lb    ci.ub 
## intrcpt          -0.0497  0.0079  -6.3202  <.0001  -0.0652  -0.0343  *** 
## IRSAD             0.0457  0.0090   5.0922  <.0001   0.0280   0.0633  *** 
## StateVictoria     0.0148  0.0126   1.1763  0.2403  -0.0099   0.0395      
## RegionRegional    0.0155  0.0125   1.2427  0.2149  -0.0090   0.0401      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
## 
## Mixed-Effects Model (k = 331; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
##  380.8468  -761.6936  -751.6936  -732.7438  -751.5067   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0028 (SE = 0.0004)
## tau (square root of estimated tau^2 value):             0.0533
## I^2 (residual heterogeneity / unaccounted variability): 63.02%
## H^2 (unaccounted variability / sampling variability):   2.70
## R^2 (amount of heterogeneity accounted for):            37.81%
## 
## Test for Residual Heterogeneity:
## QE(df = 327) = 876.0617, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 327) = 34.7422, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt           0.0372  0.0060  6.1879  <.0001  0.0254  0.0490  *** 
## IRSAD             0.0592  0.0070  8.5049  <.0001  0.0455  0.0729  *** 
## StateVictoria     0.0505  0.0101  5.0069  <.0001  0.0306  0.0703  *** 
## RegionRegional    0.0256  0.0096  2.6728  0.0079  0.0068  0.0445   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
res$Variable[c(2:3,5:6,8:9,11:12)]<-""


print_table(res, align = c("l", "l",rep("c",5)),  colwidth = c(2,3,rep(4.4,5)))
## <caption> (\#tab:)</caption>